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Presented is a semi-analytic model of magnetized liner inertial fusion (MagLIF). This model accounts for 
several key aspects of MagLIF, including: (1) preheat of the fuel (optionally via laser absorption); (2) pulsed- 
power-driven liner implosion; (3) liner compressibility with an analytic equation of state, artificial viscosity, 
internal magnetic pressure, and ohmic heating; (4) adiabatic compression and heating of the fuel; (5) radia¬ 
tive losses and fuel opacity; (6) magnetic flux compression with Nernst thermoelectric losses; (7) magnetized 
electron and ion thermal conduction losses; (8) end losses; (9) enhanced losses due to prescribed dopant con¬ 
centrations and contaminant mix; (10) deuterium-deuterium and deuterium-tritium primary fusion reactions 
for arbitrary deuterium to tritium fuel ratios; and (11) magnetized a-particle fuel heating. We show that 
this simplified model, with its transparent and accessible physics, can be used to reproduce the general ID 
behavior presented throughout the original MagLIF paper [S. A. Slutz et al., Phys. Plasmas 17 , 056303 
(2010)]. We also discuss some important physics insights gained as a result of developing this model, such as 
the dependence of radiative loss rates on the radial fraction of the fuel that is preheated. 
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I. INTRODUCTION 

The Magnetized Liner Inertial Fusion (MagLIF) 
concep^i^ is presently being investigated 
experimentally^Ti^S using the Z facility^iiS^ at San¬ 
dia National Laboratories. MagLIF is part of a broader 
class of concepts referred to collectively as magneto- 
inertial fusion (MIF)<^“— These concepts seek to 
significantly reduce the implosion velocity and pressure 
requirements of traditional inertial confinement fusion 
(ICF)^“— by using a magnetic field to thermally 
insulate the hot fuel^ from a cold pusher and to increase 
fusion product confinement. 

The MagLIF concept at Sandia uses the electromag¬ 
netic pulse supplied by the Z accelerator to radially im¬ 
plode an initially solid cylindrical metal tube (liner) filled 
with preheated and premagnetized fusion fuel (deuterium 
or deuterium-tritium). The implosion is a result of the 
fast z-pinch process, where a large gradient in the ap¬ 
plied magnetic field pressure operates near the liner’s 
outer surfacei^*^ One- and two-dimensional simulations 
of MagLIF using the LASNEX radiation magnetohydro¬ 
dynamics code^ predict that if sufficient liner integrity 
can be maintained throughout the implosion, then sig- 
nihcant fusion yield (>100 kJ) can be attained on the Z 
accelerator when deuterium-tritium fuel is used and the 
accelerator’s Marx generators are charged to 95 kV to 
obtain a peak drive current of about 27 MAiii^ 

To maintain liner integrity throughout the implo¬ 
sion, liners with thick walls have been proposedi*^ and 
used experimentally.— Thick walls mitigate 
the deleterious effects of the Magneto-Rayleigh-Taylor 
(MRT) instability^-.^>adad5d^>^>^-.^ which first devel¬ 
ops on the liner’s outer surface, and then works its way 
inward, toward the liner’s inner/fuel-confining surface. 


throughout the implosion. The parameter that is typ¬ 
ically used to describe the robustness of a liner to the 
MRT instability is the liner’s initial aspect ratio 

Ao = (1) 

OlO 

where rjo is the liner’s initial outer radius and Siq is the 
liner’s initial wall thickness. Lower A^o liners are more 
robust to the MRT instability, but their use results in 
slower implosion velocities; thus, there is a tradeoff be¬ 
tween liner robustness and implosion efficiency. 

Two-dimensional LASNEX simulations predict that 
liners with Aro < 10 should be robust enough to keep the 
MRT instability from overly disrupting the fusion burn 
at stagnation— These preliminary LASNEX simulations 
predict a broad optimum in the fusion yield surrounding 
Aro ~ 6 (see Fig. 10 in Ref. [l|). Experiment a^d° on the Z 
accelerator with no initial axial magnetic field have found 
the implosion dynamics of Aro = 6 Be liners to be in good 
agreement with the 2D LASNEX predictions found in the 
original MagLIF paper— (e.g., compare Fig. 2 in Ref. M 
with Fig. 9 in Ref. [ll). By contrast, Aro = 6 Be liner 
experiments on the Z accelerator that included an initial 
axial magnetic field of about 10 T have revealed the pres¬ 
ence of a helical instability structure that is fundamen¬ 
tally 3D in naturCfi^d^ and thus could not be captured by 
the 2D LASNEX simulations; regardless, the inclusion of 
the initial axial field seems only to have improved overall 
implosion stability—^ Moreover, the first fully-integrated 
experimental tests of MagLIF have shown that a fuel im¬ 
plosion driven by an A^o = 6 Be liner, combined with pre¬ 
heat and premagnetization, does indeed result in fusion 
relevant conditions on the Z accelerator.— Therefore, 
we assume that the qualitative ID versus 2D behavior 
shown in Fig. 10 of Ref. [l| holds fairly true for A^o < 6, 
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and thus that MagLIF can be modeled reasonably well 
in ID for Aro ^ 6. 

The leading hypothesis for what seeds the MRT de¬ 
velopment is the electrothermal instability (ETI)ji2r2i 
Recent experiments have shown that ETI growth can be 
tamped by applying a thick dielectric coating to the outer 
surface of the metal liner, thereby delaying the onset of 
nonlinear MRT growth.— Thus, it is hopeful that the ap¬ 
plication of dielectric coatings may enable MagLIF’s 2D 
and/or 3D « ID performance to be extended to liners 
with aspect ratios significantly greater than 6. Note that 
these assumptions do not account for possible azimuthal 
variations; the effects of azimuthal variations on fuel con¬ 
finement are beyond the scope of this paper, but they are 
presently being investigated in other studiesPre¬ 
sumably, the deleterious effects of azimuthal variations 
are also mitigated by the use of low Aro liners, and per¬ 
haps by dielectric coatings as well. 

Since MagLIF calls for low Aro liners to mitigate the 
MRT instability, the resulting liner implosions are too 
slow to shock heat the fuel (i.e., the fuel heating due to 
compression is approximately adiabatic). Thus, to heat 
the fuel to fusion-relevant temperatures (>1 keV) while 
simultaneously lowering fuel convergence requirements, 
the MagLIF concept calls for preheating the fuel just 
prior to the implosion. The parameter typically used 
to describe fuel convergence is the convergence ratio 

where r^o and rg{t) are the initial and time-dependent 
radii of the fuel-liner interface, respectively. The con¬ 
vergence ratio is often evaluated at particular times of 
interest, such as the time of peak fusion power, Crp, and 
the time of minimum r^, i.e., the time of “bounce”, Crb- 
In many cases, Crp and Crb have fairly similar values, 
but differences between these two parameters can indi¬ 
cate important phenomena such as fuel collapse or fuel 
ignition. To reduce Crp and Crb requirements to < 30, 
the MagLIF concept calls for preheat temperatures of 
> 100 eV. 

At the Z facility, fuel preheating has been accomplished 
using the Z beamlet laser (ZBL)^i^ which was origi¬ 
nally coupled to the Z pulsed-power accelerator for diag¬ 
nostic purposes (e.g., radiograph y . In the first fully- 
integrated MagLIF experiments)^ ZBL provided about 
2 kJ of 532-nm light in a 2-ns pulse. Recently, ZBL has 
been upgraded to deliver about 4 kJ in a 4-ns pulse, and 
plans are in place to further increase ZBL’s delivered 
pulse energy to 6 kJ and beyond over the next several 
years. 

To keep the fuel hot during the relatively slow implo¬ 
sion of a MagLIF liner, the concept requires premagneti¬ 
zation of the fuel using an axially-aligned field. About 
10-50 T are to be supplied by external coils. This 
initial seed field is to be amplified by a factor of 10^- 
10^ within the fuel-filled volume of the imploding liner 
via magnetic flux compression. This large axial field is 


required to mitigate energy loss from the fuel due to elec¬ 
tron and ion thermal conduction. Additionally, the axial 
field should enhance a-particle confinement and heating 
of the fuel, and thus increase the overall fusion yield. 

In the first fully-integrated MagLIF experiments)^ as 
well as in liner implosion dynamics experiments,— ini¬ 
tial seed fields of up to 10 T have been supplied by a 
newly developed applied B on Z (ABZ) subsystem at 
the Z facility— Analysis of the neutron diagnostics data 
collected during the first fully-integrated MagLIF exper¬ 
iments indicates that significant flux compression did in¬ 
deed occur— Also, the ABZ capabilities at the Z facil¬ 
ity have recently been upgraded to 15 T, and plans are in 
place to further increase these fields to 30 T in the next 
couple of years— 

To elucidate some of the key physics issues relevant to 
MagLIF, we have developed a new semi-analytic model of 
the concept. This model is formulated as a system of or¬ 
dinary differential equations (ODEs) that are straightfor¬ 
ward to solve, particularly with standard software tools 
such as MATLAB®, IDL®, Mathematica®, etc. This 
model accounts for: (1) preheat of the fuel (optionally via 
laser absorption); (2) pulsed-power-driven liner implo¬ 
sion; (3) liner compressibility with an analytic equation 
of state, artificial viscosity, internal magnetic pressure, 
and ohmic heating; (4) adiabatic compression and heat¬ 
ing of the fuel; (5) radiative losses and fuel opacity; (6) 
magnetic flux compression with Nernst thermoelectric 
losses; (7) magnetized electron and ion thermal conduc¬ 
tion losses; (8) end losses; (9) enhanced losses due to pre¬ 
scribed dopant concentrations and contaminant mix; (10) 
deuterium-deuterium and deuterium-tritium primary fu¬ 
sion reactions for arbitrary deuterium to tritium fuel ra¬ 
tios; and (II) magnetized a-particle heating. This model 
has been implemented in a code called SAMM (Semi- 
Analytic MagLIF Model). Simulations using SAMM typ¬ 
ically take about 30 seconds to run on a laptop using the 
ode23 solver in MATLAB®. Using the parallel com¬ 
puting cluster at Sandia, parameter scans of about 2000 
simulations can be completed in as little as 10 minutes. 

Semi-analytic models have proven themselves useful in 
the past— Their fast run times allow the system pa¬ 
rameter space to be explored rapidly. Also, the physics 
included in semi-analytic models are often transparent 
and thus accessible by anyone, which is not always the 
case for more sophisticated codes (e.g., the LASNEX 
code^ that the original MagLIF paper— was based on). 
Finally, one of the most important aspects of develop¬ 
ing simplified models such as SAMM, is that it forces 
those involved to verify and understand the results gen¬ 
erated by more sophisticated codes. For example, when 
differences are observed between a simple model and a 
more advanced code, can those differences be explained 
in terms of the simplifying assumptions made in the re¬ 
duced model? This type of questioning often leads to new 
physics insights that would otherwise be overlooked. For 
example, during the development of SAMM, we found an 
important relationship between radiative loss rates and 
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the initial radial fraction of the fuel that is preheated; this 
finding and others are discussed further in this paper. 

The remainder of this paper is organized as follows. In 
Sec. m we present our semi-analytic model of MagLIF. 
In Sec. uni we verify this model relative to the ID results 
presented in the original MagLIF paper (i.e., Ref. [U). 
In Sec. lYl we summarize this work. Unless otherwise 
specified, all units are SI. 


II. A SEMI-ANALYTIC MODEL OF MAGLIF 


A. Overview of the model and the radial distribution of 
the driving azimuthal magnetic field 


An overview of the semi-analytic MagLIF model is pro¬ 
vided in Fig. [TJ The liner implosion, and thus the fuel 
implosion, is driven by the pressure associated with the 
azimuthal magnetic field, which is supplied by the pulsed- 
power driver (e.g., the Z accelerator). Because of the 
cylindrical symmetry, the azimuthal field in the vacuum 
region is 


B0v{r) 


2'Kr ’ 


( 3 ) 


where /xq = 47r x 10“^ H/m is the permeability of free 
space and Ii is the liner current. We assume that Bg is 
partially diffused into the liner wall, and that its distri¬ 
bution in the liner region can be described by 


Bei{r) 


fJ-oh f r-rg X 

2-Kri \n-rg) 


( 4 ) 


where the constant power j3 is found by forcing the am¬ 
plitude of Bgi{r) to drop by one e-folding within one skin 
depth, Sskin, of the liner’s outer surface; that is, we set 

Bgi{ri - Sskin) = ■ -, (5) 

Inri e 


where e = 2.71828 is the base of the natural logarithm, 
and solve for /?, which gives 


/3 


ln(l/e) 


In 


V ^‘0 ) 


( 6 ) 


The skin depth is given by 


_ / ipeTr 

Oskin — \ •) {• ) 

V 

where pe is the initial electrical resistivity of the liner 
material, and is the rise time of the driving Bg pulse 
(e.g., ^ 130 ns for the Z accelerator). Some parameters of 
interest are provided in Table HI The purpose of assuming 
a simple power-law dependence for Bgi(r), rather than an 
exponential dependence, is that it allows us to integrate 
both the total magnetic flux and the total magnetic field 
energy analytically, which enables us to obtain simple 
analytic expressions for the inductance, circuit response, 
and ohmic dissipation due to Bgi. 



FIG. 1. Schematic overview of the semi-analytic MagLIF 
model. There are three primary regions: a fuel region, a 
liner region, and a vacuum region. The system height is h; 
the thermally-insulating axial magnetic field, which is initially 
distributed uniformly over all regions, is Bz ; the radius of the 
fuel-liner interface is rg\ the liner’s outer radius is rp the re¬ 
turn current radius is r^c; the azimuthal magnetic field, which 
drives the cylindrical implosion, is Bg. Normalized profiles are 
shown for Bg in the vacuum region, Bg^ (magenta), and in 
the liner region, Bgi (orange); their analytic expressions are 
given by Eqs.[3]and[3]in the text {Bg is assumed to be zero in 
the fuel region). The liner region is further divided into mul¬ 
tiple concentric liner shells; this discretization is necessary to 
avoid overdriving the fuel. Within the fuel region, normal¬ 
ized profiles are shown for the gas pressure, Pg (blue), the gas 
temperature, Tg (red), the gas density, pg (green), and the 
radiation temperature, Tr (cyan). The pressure profile is flat 
throughout the fuel (i.e., we have made an isobaric assump¬ 
tion due to the subsonic nature of MagLIF implosions). The 
gas temperature and density profiles thus have an inverse de¬ 
pendence to one another; their analytic expressions are given 
by Eqs. 1105111091 in the text. The radiation temperature is 
nearly constant across the fuel region. The fuel region is fur¬ 
ther divided into a hot spot region from r = 0 to and a 
cold dense shelf region from rg to Tg. The gas temperature 
in the shelf region is equal to the radiation temperature (i.e., 
the fuel material and the radiation field are in thermodynamic 
equilibrium in the shelf region). The shelf region erodes away 
throughout the implosion, until rg = Tg, due to thermal trans¬ 
port from the hot spot to the shelf. The shelf region is only 
present if the fuel is preheated from r = 0 to r < r^. 


TABLE 1. Parameters for Bgi{r) (where Tr = 130 ns). 


Material 

Pe [nD-m] 

^skin [mIh] 

P 

Li 

92.8 

110.6 

3.683 

Be 

36 

68.9 

6.239 

A1 

28.2 

60.9 

7.118 


B. Circuit model for the Z pulsed-power driver 

The liner current in Eqs. [3] and |4] can be specified to 
drive the simulation directly, or it can be derived from 
a circuit model. For the circuit model option, we use 
the circuit illustrated in Fig. [21 which has been shown 
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4 L 


Lq 



FIG. 2. An equivalent circuit model for the Z pulsed-power 
accelerator, where ipoc{t) is the open-circuit voltage used to 
drive the simulation (see Fig. [3]), Zo=0.18 fl, L=8.34 nH, 
(7=8.41 nF, and Lo«5 nH. The elements Lv{t) and Lic{t) are 
given by Eqs. [16] and in the text and they represent the 
coupling of the circuit to the dynamic volume of the model 
illustrated in Fig. [T] 


and thus 


L 


^oc 

z 


The liner current is 


(9) 


Z - Z CiPc ^c/Rlossi 


and thus 


Zc 


Z Z ^c/ Rloss 

c 


( 10 ) 


( 11 ) 


The voltage required to supply Bg flux to the model’s 
volume illustrated in Fig. [1] is 




( 12 ) 



FIG. 3. Example open-circuit voltage waveform for the Z 
pulsed-power accelerator. This waveform was constructed us¬ 
ing twice the forward-going voltage measured at the Z accel¬ 
erator’s vacuum-insulator stack. 


Since the axial current density in the liner wall goes to 
zero at r = rg, which is implied by our assumption for 
Bgi{r), then by the integral form of Faraday’s law 

j(E-dl = -4>e, (13) 

where the path integral curve c is around the model’s 
volume illustrated in Fig. jT] E is the electric field vector, 
dl is an infinitesimal path element vector along curve 
c, and is the total Bg flux in the model’s volume 
illustrated in Fig. |l] we know that we can write iprc also 
as 

^rc = ^e = ^ev + ^eu (14) 

where and ^gi are the total Bg fluxes in the vacuum 
and liner regions, respectively. For the vacuum region, 
we have 

— l^vZ 4“ Z^Z? (1^) 

where is the standard coaxial vacuum inductance^ 


to adequately represent the Z accelerator coupled to z- 
pinch loadsi^ With this model, simulations are driven 
by an open-circuit voltage, <~Poc(t)- This voltage is twice 
the forward-going voltage at the vacuum-insulator stack 
on Z, and it can be constructed from the data of previous 
Z experiments using Eq. 4 in Ref. Ib^. An example (fiocit) 
waveform for Z is shown in Fig. |3| 

This circuit model also requires Riossit) to be speci¬ 
fied. This is a resistance used to model shunt current 
losses on Z, and it can be derived from previous Z exp er- 
iment data using the methods discussed in Refs. l62Hd4L 
For simplicity, we will take Riossit) to be large and con¬ 
stant (>10 11) throughout this paper, which essentially 
eliminates current loss from the results presented herein. 

To solve the circuit shown in Fig. [21 we begin by writing 
the voltage across the capacitor 


— ^oc -^oZ RRf 


( 8 ) 


and thus 


Ly{t) = ^In 


Mt)\ ’ 


Lvit) 


MqZ / rA 

27r \ri) 


(16) 


(17) 


For the liner region, we have 

fJ-ohliin - rg) 


^ei = 
and thus 
^ei = 


rri 

h / Bgi{r) ■ dr = 

Jr„ 


27rr/(/3-I-1) 


(18) 


/ioh 


27r(/3-bl) _ 


z 1- 


ri 


TgB rg 

rf 


n 


(19) 


It is convenient to define an effective inductance that can 
be used to describe the circuit response due to the Bg flux 
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in the liner region. From Eq. 1191 we define this effective 
inductance as 


and thus 


Licit) = 


Licit) = 


fioh 


27r(/3 + l) 
fJ-oh 


1-^ 

n 


Tgn 


( 20 ) 


( 21 ) 


27r(/3 + 1) V n , 

With these expressions, Eq. [TOlcan be rewritten as 

= Licii + Lich- (22) 


Finally, combining Eqs. [T2l El El and [22l and solving 
for Ii gives 


ii 


‘fc ~ hiLv + Lie) 

Lq + Ly + Lie 


(23) 


C. Drive energetics and ohmic liner heating 


From Eq. El the power going into just the Bg field is 
F’bo = Lbd = iLy + Li)iili + -^iLv + Li)lf. (28) 

To obtain an expression for Pkin , we consider the case of 
a perfectly conducting liner, where Pq, = 0, and Lie = 
Li = 0 because the liner current is a surface current at 
r = ri (i.e., /3 —>■ oo)i^ This gives 

PkiniP oo) = PeM — PBb = -Lylf, (29) 

which is an upper bound for Pkin- From numerical tests 
with finite conductivity, we find that to within a few 
percent over most of the implosion 


Pkin ~ PkiniP -)■ oo) = ]^Lyl‘f. (30) 

(Below, in Eq. 1481 we provide a more precise definition 
and explanation of Pkin-) Finally, the power associated 
with ohmic dissipation is found from 


It is important to realize that the definition for Lie 
above is not the same as the definition for standard in¬ 
ductance; the standard inductance is found from the total 
magnetic field energy in the liner region^ 


Pq = \PeM — Pbb — Pkm \ - (31) 

Substituting Eqs. El HH and El into Eq. El gives the 
approximate (and lower bound) solution 


1 

Pbb, = ^Lilf = hj ^ • 27rr • dr, 
which gives 

j- ^ A^ofe(2/3r; +ri+ rg)(r; - Vg) 

* 47rrf (/3-I-l)(2/3-I-1) 

and thus 

f _ ^^ohirg + Pri)irgri - rgj-i) 

27rrf(/3 + l)(2^ + l) ' 


(24) 


(25) 


(26) 


Note that with this standard inductance, $ 61 ; ^ Lili + 
Lili. Differences exist between Lie and Li because the 
liner current is distributed radially across the liner region. 
The distributed current is a result of our assumed Bgiir), 
which we chose in an attempt to model the magnetic 
diffusion and ohmic dissipation known to occur in thick- 
walled liner implosions 

These two definitions for inductance provide us with 
a clean and simple way to express the following powers 
being delivered to the volume of Fig. [H (a) the total 
electromagnetic power, Pem', (b) the power going into 
just the azimuthal magnetic field, Pbb] (c) the power 
going into the kinetic motion and compression of the fuel, 
liner, and axial magnetic field, Pkin] and (d) the ohmic 
heating power, Pq. 

The total electromagnetic power is the sum Pem = 
Pbb + Pkin + Pq] it can also be expressed in terms of the 
drive circuit as 


Pem — ^eli — iLy -L Lic)iili -L iLy -I- Lic)!^ . (27) 


Pq 



Li)iili + 



(32) 


Note that Eqs. El EU and [32] are valid only if the 
pulsed-power generator can supply azimuthal flux to the 
liner wall faster than the flux can be dissipated ohmi- 
cally, and thus Sskin < ri—Vg. In typical MagLIF experi¬ 
ments, this condition is essentially always met due to the 
use of electrically-conductive, thick-walled liners, and be¬ 
cause the pulsed liner current changes rapidly with time 
throughout almost the entire experiment.— In numer¬ 
ical tests, we have found that Eqs. E and E] generate 
ohmic dissipation rates that agree well with those cal¬ 
culated by full radiation magnetohydrodynamics simula¬ 
tions. 

Accounting for Pq is important because without it, 
the liner region remains cooler and more compressible 
than expected throughout the implosion (i.e., on a “lower 
adiabat”). Upon stagnation, an overly compressed liner 
region will drive the fuel to a pressure that is too high, 
and thus result in an overly optimistic fusion yield. 


D. Liner dynamics and compression 

To account for liner compressibility, we first divide the 
liner region into Nis > 20 concentric thin liner shells and 
write an equation of motion for the interface between 
each shell (internal interfaces), as well as for the fuel-liner 
interface and for the liner-vacuum interface (external in¬ 
terfaces) . To each internal interface, we assign a mass of 
mis = TBi/Nis, where mi is the total liner mass; to each 
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external interface, we assign a mass of misl2. There are 
Nii = Nis + 1 liner interfaces in total. The radial po¬ 
sitions of the liner interfaces, ri^i, are distributed from 
rg = ri^i—i to r; = The initial radial position of 

interface * -I- 1 is given by 

n,i+i(io) = \jrliih) + mis/{Trhpio) 

(* = l,2,...,iVH-2), (33) 

where pio is the initial mass density of the liner material. 
The center of mass radial position for liner shell s is given 

by 


n.sit) = 


rli=sit) + /2 


(s=l,2,...,iVi,). (34) 


Since there are Nn interfaces, there are Nn equations of 
motion. The equations of motion for the internal liner 
interfaces are given by 


Pl.s — i — l Pl^S—i r, 1 

ri,i = -• 27rr;_i • h 


mis 


(i = 2,3,...,7 Vh-1), (35) 


where pps is the effective pressure in liner shell s. Each 
pps is comprised of material pressure, Pmi,s, azimuthal 
magnetic field pressure, pbbi.s, axial magnetic field pres¬ 
sure, Pb^i,s, and pseudo-pressure due to artificial viscos- 
ityi ®,s- The equation of motion for the fuel-liner inter¬ 
face is 


Pa+PS,, -Pi,s=i 

mis 12 


■ 2'!Trg ■ h, 


(36) 


where pg is the gas pressure in the fuel region and pg^g is 
the pressure due to the average axial magnetic field in the 
fuel region. The equation of motion for the liner-vacuum 
interface is 


h 


Pl,s=Ni, - PBei^ - Pb,^ 
mis (2 


■ 2'Kri ■ h. 


(37) 


where PBgig is the pressure of the azimuthal magnetic 
field evaluated at the liner-vacuum interface {r = ri) and 
Pb^^ is the pressure due to the average axial magnetic 
field in the vacuum region. 

The magnetic field pressures are simply 

PBe-R.{r) = (38) 

where C is the component (either 9 or z) and TZ is the 
region (either g for the fuel/gas region, I for the liner 
region, or v for the vacuum region). We also note the 
following. Since we have made an isobaric assumption 
in the fuel, and since, for MagLIF, the gas pressure typ¬ 
ically dominates the magnetic field pressure in the fuel, 
we simply use the average axial magnetic field in the fuel 


B 


zg ~ 


9 ’ 

n I g 


(39) 


for Pg , where ^zg is the total axial magnetic flux in 
the fuel. By contrast, in the liner region, we assume 
Bziir)/pi(r) = const., and thus for PB„i,s, we use 


Bzi,siri,s{t)) = Bzi^ 

pi 


where 


Bzi = 


4) 


zl 


(40) 


(41) 


TT{rf - r-2) 

is the average axial magnetic field in the liner region, ^zi 
is the total axial magnetic flux in the liner region 

mi 


n{rf - r^g)h 

is the average mass density in the liner region 

mis 


Pl,s — 


Vl. 


(42) 


(43) 


is the mass density of shell s, and 

Vps = Tl'{flA=s+l ~ ^l,i=s)^ ( 44 ) 

is the volume of shell s. For pBgi , we evaluate our analytic 
expression for Bgi (r) (see Eq. |4]) at the center of mass of 
each liner shell, n,s(i) (see Eq.jM]). For we use the 
average axial magnetic field in the vacuum 


, 2^" 2V (45) 

[rrc - rf) 

where ^zv is the total axial magnetic flux in the vacuum 
region of Fig. [TJ Finally, for PBgi^, we use the azimuthal 
magnetic held evaluated at r = r/ 

Poll 


Bgiv = Bgi(ri) = Bgy{ri) = 


2'Kri 


(46) 


and thus the resulting azimuthal held pressure at r=ri is 

Polf 


PBgi„ = 


BL 


2p,Q 


8k\^ 


2 ■ 


(47) 


With the magnetic held pressures and liner interface 
equations of motion thus dehned, we can revisit Eqs. 127b 
1321 Specihcally, we can now obtain a more precise ex¬ 
pression for the kinetic power 


Pkin = ^ S(j)Bg)t ■ ri,i ■ ‘^KTp^ ■ h, 


(48) 


where S{pBg)i is dehned as the difference in azimuthal 
magnetic held pressure on either side of interface i. This 
expression describes the rate at which the driving Bg held 
does work on each interface, i, and thus the rate at which 
the Bg held does work on the fuel and liner materials (via 
material compression and acceleration), as well as on the 
axial magnetic held, Bz (via hux compression). This ex¬ 
pression for Pfem can be used to replace the approximate 
expression given in Eq. 1301 and then from Eq. 1311 we can 
obtain a more exact expression for Pq, thus replacing 
Eq. 1321 However, in practice, both Eos. [501 and H51 work 
hne, since typically, they are within a few percent of each 
other over most of the implosion. 


fJournal Reference: Phys. Plasmas 22 , 052708 (2015); http://dx.doi.Org/10.1063/l.4918953 
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TABLE II. Nuclear charge, atomic mass number, and cold 
curve fitting parameters for Eq. [50]for various liner materials. 


Material Znuc 

A 

Pio [kg/m’^j 

Al 

A2 

71 

72 

Li 

3 

6.94 

534 

11 X 10'^ 

3.999 

1.9 

1.18 

Be 

4 

9.012 

1845 

130 X 10® 

3.9993 

1.85 

1.18 

Al 

13 

26.98 

2700 

76 X 10® 

3.9 

7/3 

5/3 


E. Liner equation of state, energetics, and ionization 


range of 1-3 works well for our purposes. Note Eq. [51] 
states that the artificial viscosity contributes to the shell 
pressure only during times of compression. The rate of 
change of the internal thermal energy of liner shell s is 
thus given by 


E 


/,s — 


+ {Pr + Pc + Pq. — Pbb) /Nls 


(53) 


To model the internal material pressure of the liner 
shells, we use the expression 


— ( 1 -L —Ehl 

Pml,s — P0\Pl,s) P Q -tr 
O VLb 


(49) 


The first term, po{pi,s)j is the zero temperature “cold 
curve” for the liner material, which is a function of only 
the mass density, pi^s- To model the cold curves of vari¬ 
ous liner materials, we use the functional form of Birch- 
Murnaghan22,“— 





72-1 


X 


1 + - [^2 



(50) 


where pio is the liner material’s mass density at zero tem¬ 
perature and zero pressure, and where Ai, A^^ 7i, and 
72 are fitting parameters. For lithium (Li), beryllium 
(Be), and aluminum (Al), we find reasonably good fits 
to SESAME equation of state datai^ using the values 
provided in Table im 

The second term on the right hand side of Eq. |49|is the 
ideal gas thermal contribution to the equation of state, 
where Ei g is the total thermal energy of the liner ma¬ 
terial in shell s. With finite liner temperature, Eps can 
change due to adiabatic compression and/or expansion. 
Additionally, Ep^ can increase because of artificial vis¬ 
cosity and ohmic heating. We have found that includ¬ 
ing artificial viscosity is useful because it dampens the 
spring-like reverberations that occur in the liner’s wall 
thickness throughout the implosion. The artificial vis- 
cosity model that we use is similar to that discussed in 
Ref. 1^. Our pseudo-viscous pressure for shell s is given 
by 


= 


alppsirpi=s+i - ri,i=s)^ for rpi=B > rpi=s+i 
^0 for fpi=s < rpi=s+i, 

(51) 

where a, is an arbitrary coefficient of 0{1) that specifies 
the length scale of the artificial viscosity in multiples of 
the radial thickness of shell s 


^l,s —s-t-1 rpi—B' 


(52) 


where 

Vps = 2TTh (n,i=s+in,i=s+i - rpi=sh.i=s) ■ (54) 

The first term on the right hand side of Eq. |S3I is due to 
adiabatic compression/expansion, while the second term 
is due to artificial viscous heating during compression 
only. The quantities Pr and Pc represent the radiative 
and thermal conduction loss powers from the fuel region, 
respectively; they will be described in more detail below. 
For simplicity, the powers Pr, Pc, and Pq are all assumed 
to be fully absorbed by the liner and evenly distributed 
throughout the liner, hence they are each divided by Nis 
in Eq. |5l2^ The liner is permitted to cool through adia¬ 
batic expansion as well as blackbody radiation. We use 
blackbody radiation because the liner is optically thick 
for typical liner materials (Li, Be, Al), wall thicknesses 
(>100 A^m): and average liner temperatures (<100 eV). 
The blackbody power is given by 

Pbb = ■ 2nri ■ h, (55) 

where a = 5.67 x 10“® W/(m^-K'^) is the Stefan- 
Boltzmann constant. For simplicity, energy losses from 
the liner due to Pbb are evenly distributed throughout 
the liner, hence Pbb is also divided by Ws in Eq. |53l 
Note that in Eq. 1551 we use the average liner tempera¬ 
ture, which is given by 

T = ^ _ 2 uAJ2sPi,s .gg, 

' 3 Nik 3Nat{l + Zi)k 3mi{l + Zi)k' ^ ’ 

where k = 1.38 x 10“^^ J/K is the Boltzmann con¬ 
stant, El is the total thermal energy of the liner, Ni = 
Nat{f + Zi) is the total number of particles in the liner 
(ions and electrons), Nat = mi/{uA) is the total number 
of atomic nuclei in the liner, u = 1.66 x 10“^^ kg is the 
unified atomic mass unit, A is the atomic mass number 
of the liner material (see Table Hill . and Zi is the average 
ionization state of the liner. We approximate the average 
ionization state using2& 

Zi = min (^O^JfpkeV, Znu^ , (57) 

where Znuc is the atomic number (nuclear charge) of the 
atoms comprising the liner material (see Table m. The 
average liner temperature in keV is simply 


The overall fusion calculations are not very sensitive to 
the exact value chosen for Oq. We find that Oq in the 


Tl,keV 


m 

Qe 


X 10 


-3 


(58) 


1 Journal Reference: Phys. Plasmas 22 , 052708 (2015); http://dx.doi.Org/10.1063/l.4918953 
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where qe = 1.6 x 10“^® C (or J/eV) is the charge of an 
electron. 

For reasonable fusion calculations, we have found it 
necessary to discretize the liner region into concentric 
thin shells, as described above, because even a compress¬ 
ible cylindrical slab-like model for the liner region (i.e., 
no internal interfaces, just the fuel-liner and liner-vacuum 
interfaces connected by the equation of state described 
above) results in very optimistic fusion yields due to very 
optimistic fuel compression. This occurs because in the 
cylindrical slab-like model, too much mass is assigned to 
the fuel-liner interface. The only way to reduce the mass 
at the fuel-liner interface is to use a finer discretization. 
We find that solutions tend to converge when approxi¬ 
mately 20 liner shells are used. 


F. Fuel ionization, energetics, and adiabatic heating 


Within the fuel region, we assume adiabatic compres¬ 
sion and a radially constant pressure profile (an isobaric 
assumption). In reality, there will be small radially de¬ 
pendent pressure waves that reverberate within the fuel; 
however, the slow liner implosion, which on average is 
subsonic relative to the fuel’s sound speed,— implies that 
our adiabatic approximation is well justified. This also 
allows us to assume equal ion and electron temperatures, 
which we do throughout this semi-analytic model. 

To model the energetics in the fuel region, we assume 
that the fuel is an ideal gas, and thus the isobaric particle 
pressure is related to the total thermal energy of the fuel, 
Ethg, by 


2 Ethg 

where 


(59) 


Vg = ‘^^fgh (60) 

is the volume of the gas. Furthermore, the sum of the 
thermal energy and the ionization energy gives the total 
energy of the fuel 


Eg — Ei 


thg 


Ei, 


(61) 


We have found that accounting for the ionization energy 
is important for MagLIF in cases where the available pre¬ 
heat energy is low (^ 100 J). For example, about 50 
J would be required to fully ionize the deuterium fuel 
used in recent experiments at the Z facilityj^ This is 
a substantial fraction of the 100-300 J of preheat en¬ 
ergy that is thought to have coupled to the fuel in these 
experimentsi^ 

For the ionization energy of particle species s in the 
fuel, we have found good fits to tabulated data— using 
the analytic expression 


Eion,s — 13.6 * q^ 




i^Zjiy^CS Zg -\- 1 ^ 


Vc 


C = 2.407, 


(62) 

(63) 


where Znuc,s is the atomic number (nuclear charge) of 
atomic species s, and Zg is the average ionization state 
of species s. Here we again make use of the approximate 
formula— 


where 


Zg = min I 20\/T, 


- g,keV 5 ^nuc,s I 5 


- _ kTg 3 

q.keV — ^ 

Qe 

rp _ 2 Ethg 

{N, + Ne)k 

N, = Y,Ng 

S 

Ne = ZgNi 


(64) 

(65) 

( 66 ) 

(67) 

( 68 ) 

(69) 


and where Tg^keV and Tg are the mean temperature of 
the entire fuel region in units of keV and K, respectively; 
Ni and W are the total number of ions and electrons 
in the fuel, respectively; Ng is the total number of ions 
of atomic species s in the fuel; and Zg is the average 
ionization state of the fuel. Finally, summing over all 
species gives 


Eiong — ^ ^ Eton,s ■Ng. (70) 

S 

The purpose of accounting for various particle species, 
s, is so that we can study the effects of prescribed levels of 
dopants and/or contaminants (“mix”) in the fuel region. 
Also note that since ionization depends on temperature 
and temperature depends on ionization, this would re¬ 
quire an iterative solution within our overall ODE system 
solve. However, since we are only roughly estimating the 
ionization, and since it changes relatively slowly, we avoid 
having to implement an iterative solution by holding the 
ionization fixed during a system solve and updating it 
only between successive system solves (which is typically 
every ^ 100 ps). 

With ionization established, we note the following sim¬ 
ple relationships, which will be used throughout the re¬ 
mainder of this manuscript: 


rUg = y^uAgNg 

S 


= -nidNd + mtNt -f- ^ uAgNg 
s^d,t 

(71) 

rhi = nig/Ni 

(72) 

Pg = mg/Vg 

(73) 

fii = Ni/Vg 

(74) 

fie = Ne/Vg 

(75) 

fig = Ng/Vg, 

(76) 


^Journal Reference: Phys. Plasmas 22 , 052708 (2015); http: //dx. doi . org/10.1063/1.4918953 
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where nig is the total mass of the fuel; As is the atomic 
mass number of species s; md = 3.34 x 10“^^ kg is the 
mass of a deuteron; m* = 5.01 x 10“^^ kg is the mass 
of a triton; Nd and Nt are the number of deuterons and 
tritons in the fuel, respectively; ifii is the average mass of 
an ion in the fuel; pg is the average mass density in the 
fuel; hi and he are the average ion and electron densities 
in the fuel, respectively; and h^ is the average ion density 
in the fuel for species s. 

The dynamics of the total fuel energy are described by 
Eg = PpdV T Pph T Pa Pr Pc Pends 


that is present in case the laser propagates completely 
through the imploding region. If the fuel reservoir were 
not present, the laser could hit and ablate high-Z mate¬ 
rial from various target and/or electrode surfaces, which 
could spray back into the fuel of the imploding region, 
providing a source of contaminant mix. 

To account for the energy absorbed in each region (i.e., 
the LEG region, the imploding region, and the beam 
dump region), we use the inverse bremsstrahlung model 
outlined in Ref. [J By assuming an ionization state for 
the fuel while it is interacting with the beam; for exam¬ 
ple, we use the fully-ionized condition given by 


where Ppdv is the adiabatic heating rate, Pph is the fuel 
preheating rate, Pc, is the heating rate due to a-particle 
energy deposition, Pr is the radiative cooling rate. Pc is 
the cooling rate due to electron and ion thermal conduc¬ 
tion, and Pends is the energy loss rate due to fuel escap¬ 
ing out of the top and/or bottom ends of the imploding 
cylinder. The adiabatic heating rate is given by 

4 

Ppdv ~ Pg^g ~ '^Pthg'^g/'^g- (7^) 


G. Fuel preheating (optionally via laser absorption) 


The fuel preheating rate can be described simply using 
a square pulse 


Pph — * 


0 

-®P'*/Tpfc 

0 


for t < tph 

for tph ^ t ^ {tph 'Tph) 
for t>{tph + Tph), 


(79) 


where Pph, Tph, and tph are the preheating energy, pulse 
length, and turn-on time, respectively. 

The experimental MagLIF program at Sandia National 
Laboratories has been investigating the approach of pre¬ 
heating the fuel using the Z beamlet laser.— In general, 
MagLIF does not require that the preheating be done 
with a laser. One could imagine finding a way to pre¬ 
heat the fuel using some of the energy supplied by the 
pulsed-power driver, which may enable more total pre¬ 
heat energy, as well as more efficient preheating of the 
fuel. Nevertheless, should the preheating be done with 
a laser, we can account for the beam propagation and 
deposition analytically. 

An analytic description of laser propagation and de¬ 
position is particularly useful for quickly evaluating the 
various MagLIF target designs presently being consid¬ 
ered for use at the Z facility. For example, MagLIF tar¬ 
gets (see Fig. 2 in Ref. [Tsh have consisted of a cylindrical 
laser entrance channel (LEG), filled will fuel, that resides 
axially between the bottom of the laser entrance win¬ 
dow (LEW) and the top of the imploding region {z = h 
in Fig. [T]). We define the axial length of this channel 
as Azlec- Furthermore, a beam dump has been used, 
residing below the imploding region (below 2 ; = 0 in 
Fig. [T|). This beam dump is simply a reservoir of fuel 


Zb 


iX 


sNs-, 


(80) 


and by ignoring thermal conduction and hydrodynamic 
motion (i.e., we assume fast isochoric laser heating), the 
heating of the fuel by laser light can be described by2^ 


deb 

dt 

K{z,t) 


dib I 

— = -K(z,t)Ib 

az 


^eA 


C UJt 




(81) 

(82) 


where eb{z,t) is the energy density in the fuel that is 
interacting with the beam, Ib{z,t) is the beam intensity 
in the fuel, k{z, t) is the absorption coefficient, Vci is the 
electron-ion collision frequency, c = 3 x 10® m/s is the 
speed of light in vacuum, Wp is the plasma frequency 
in the fuel, and LOb is the laser frequency. The plasma 
frequency is given by2S 


UJp — 



(83) 


where rric = 9.1 x 10“®^ kg is the mass of an electron, 
and eo = l/(/roc^) is the permittivity of free space. The 
laser frequency is given by 


Wfe = 27rc/Ab, 


(84) 


where \b is the laser wavelength. For preheating the fuel 
with the Z beamlet laser— at the Z facility, \b = 532 nm. 
The electron-ion collision frequency is given by^ 

_ 472//!]^ (^s^Lc.JlelnA 

^ei - n _ 3/„ 1 

{4^eo} 3v^(fcr,) 

where we have accounted for multiple ions species, s; Tb 
is the temperature of the fuel that is interacting with 
the beam; and In A is the Goulomb logarithm (see Ap¬ 
pendix |^. Next, using 


2 eb{z,t) 

3 TLi + fie’ 


( 86 ) 


we factor out the eb{z,t) dependence from Eq. [51] and 
evaluate the remaining factors at t = tph', we define the 


f Journal Reference: Phys. Plasmas 22 , 052708 (2015); http: //dx. doi . org/10.1063/1.4918953 
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result as 






(87) 


This implies that the only e{z, t) dependence in n{z, t) is 
that shown explicitly by substituting Eq. |86] into Eq. [85] 
and then Eq. |85| into Eq. |82l However, there is also an 
implicit e{z,t) dependence in the calculation of In A in 
Eq. |85l which we neglect for simplicity. For the kT de¬ 
pendence in In A, we simply use the average energy per 
particle required to meet our full ionization assumption. 
That is, we use Eq. [51] with Zs = Znuc.s, and thus we 
use kT = 13.6 • ge • Es in the In A calculation of 

Eq.lia 

The purpose of factoring out £b{z,t) in Eqs. ETHEZl is 
to make it clear that, for this subroutine of calculating 
laser absorption and propagation, we are simply freezing 
all dynamic variables [other than £b{z,t) and J{,(z,t)] to 
their values at the time when the laser is first applied (i.e., 
our assumption of fast isochoric laser heating). With this 
done, Eg. [5T] becomes 


d£b 

dt 


= -k{tph) ■ £b • Ibiz,t), 


which has an exact solution given byi 


eb{z,t) = £b{zLEW,t) 1 
Ib{z, t) = Ib{zLEW, t) 1 


ZlEW — Z 

zlew — zf{t) 

ZLEW — Z 


-,2/3 


Zlew — Zf(t)_ 


2/3 


where 


( 88 ) 


(89) 

(90) 


£b{zLEWjt) = 


k(tph) ■ Ib{zLEW, t) ■ [t — tph) 


-,2/5 


(91) 

= ZLEW - o ■ i • (t - tph), (92) 

^ ^b\EW 5 ^/ 


and where zlew = h+AzLEC is the axial location of the 
LEW (where the beam first interacts with the fuel), and 
Zf{t) is the axial location of the laser heating front as the 
beam “bleaches” its way through the fuel. Note that we 
have written this solution for a downward propagating 
beam, i.e., in the —z direction, and that this solution is 
applicable only for tph <t< (tph + Tph) and Zf(t) < z < 
Zlew- 

The energy deposition rates in the LEG region, in the 
imploding/fusion region, and in the beam dump region, 
are given by 


Plec = ttI ■ [htyZLEW, t) — Ib{h, t)] (93) 

Pph = T-rf ■ [Ib{h,t) - /h(0,t)] (94) 

Pdump = T-rl ■ Ibi0,t), (95) 


respectively, where Vb is the radius of the beam and 
Ib{zLEW,t) = Pin/{Trl). Here, is the laser power 
that makes it through the LEW and enters the fuel. 


For fusion calculations, we only need the energy de¬ 
posited in the imploding region from z = 0 to h, hence 
we only need Pph as given by Eq. |94| (or by Eq. |79| if 
ignoring laser absorption). We assume that the energy 
deposited in the imploding region is instantly distributed 
uniformly in both the axial and radial directions. In the 
axial direction, this is a reasonable assumption since ther¬ 
mal conduction is uninhibited in this direction (i.e., the 
applied magnetic field that undergoes flux compression, 
Bzg, is axially aligned, and thus inhibits thermal trans¬ 
port only in the radial direction). Furthermore, fully- 
integrated 2D radiation magnetohydrodynamics simula¬ 
tions have shown that the preheat energy redistributes 
axially on a timescale that is fast compared to the implo¬ 
sion time In the radial direction, full radiation magne¬ 
tohydrodynamics simulations show that a pressure wave 
(a “blast wave”) is generated by the rapid preheating of 
the fuel, which expands radially outward from approxi¬ 
mately the beam radius until it hits and reflects off of the 
liner’s inner surfaceJ^ Subsequent pressure waves then 
reverberate within the fuel, but settle to a reasonably 
isobaric state on a timescale that is short relative to the 
overall implosion time. 

Note that even though the pressure (overall energy 
density) in the fuel is assumed to be distributed uniformly 
in both the axial and radial directions, the fuel tempera¬ 
ture and density profiles are assumed to be uniform only 
in the axial direction. Moreover, the radial temperature 
and density profiles are significantly affected by the ra¬ 
dial extent of the fuel that is initially preheated relative 
to the radial extent of the overall fuel region. This is 
discussed in more detail below. 

For the radius of the beam, Vb, we ignore laser-plasma 
interactions and use the beam radius at the axial mid¬ 
point of the imploding region. Furthermore, we assume 
that the beam follows the focusing and defocusing cones 
described by^ 


"^ = 2 


Azb 

L77# 


1 -I- CbAzb 


(96) 


where Azb is the distance from the beam’s focal plane 
to the axial position of interest, //# is the beam’s “/- 
number”, 0 bf is the diameter of the beam in the beam’s 
plane of best focus, and Cb = 428.6 m“^ is a fitting pa¬ 
rameter. Note that the second term in Eq. [Ml is sig¬ 
nificant only in regions close to the focal plane; it is a 
first-order corrective term for describing the “waist” of 
the beam. As an example of using Eq. 1961 we consider 
the experiments of Ref. [3 where the //lO Z beamlet 
laser was focused to a 250-^m spot size, roughly 3.5 mm 
above the LEW; thus, with Azb = 3.5 mm. Eg. (Ml gives 
a spot size on the LEW of 450 fim, as quoted in Ref. [T^ . 
Additionally, the length of the LEG was about 2.1 mm, 
and the length of the imploding region was about 7.5 
mm; thus, with Azb = (3.5-|-2.1-|-7.5/2) = 9.35 mm. 
Eg. [Ml gives rj, = 492 /rm. For reference, VgQ was 2.325 
mm. 

Note that we have not accounted for absorption or 


1 Journal Reference: Phys. Plasmas 22 , 052708 (2015); http://dx.doi.Org/10.1063/l.4918953 
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scattering due to the laser entrance window. For this, 
we rely on separate laser-only experiments that mea¬ 
sured transmission through foils similar to those used for 
LEWs on MagLIF experiments; this was done for vari¬ 
ous spot sizes on the foil and will be presented in a future 
publicationj^i^ For now, we simply note that 

Pin = TlEW ■ Plaser, (97) 

where Tlew is the transmission through the LEW and 
Plaser IS the laser power on the vacuum side of the LEW. 
We comment that modeling, simulating, and experimen¬ 
tally diagnosing laser transmission, propagation, and ab¬ 
sorption in such a system is a very challenging problem, 
even for much more sophisticated simulation codes; these 
are presently very active areas of research within the 
MagLIF program and elsewhere. 


H. Fuel heating via magnetized a-particle energy 
deposition 

The next term to calculate in Eq. [77] is the heating 
rate due to a-particle energy deposition, P^. To do so, 
we note that the a particles are produced only by DT 
reactions and that the energy carried by each a particle 
is 


Qa = 3.5 X 10^ • (?e. (98) 

Following Ref. [H, the fraction of Qa that is deposited in 
our magnetized cylindrical fuel is calculated using 


/« 


Xg + X\ 

1 + 13a;a/9 + x\ ’ 


where 


/g 


b 




+ ) 

3 \la V962 -p 1000/ 

2 3 maVao{kTgf^^ 

{47reo}-- 

4v^ fieZlqtmJ^lnK 


TaL 

ZaQe ^zg 


(99) 

( 100 ) 

( 101 ) 

( 102 ) 

(103) 


and where Zg is the mean free path of an a parti¬ 
cle, TOg = 6.64 X 10“^^ is the mass of an a particle, 
Ugo = \/‘^Qalma is the birth velocity of an a particle, 
Zg = 2 is the charge of an a particle, VaL is the Lar- 
mor radius of an a particle, and the In A calculation is 
described in Appendix]^ The a particle heating rate is 
then given by 


Fg = iVdtQg/g, (104) 


where is the primary DT reaction rate, which is given 
below (see Eg. 11591) . 


I. Model of fuel hot spot and dense outer shelf regions 


For the remaining terms in Eq. 1771 which all describe 
energy losses from the fuel, and for describing magnetic 
flux loss due to the Nernst thermoelectric effect, we need 
to assume something about the temperature and den¬ 
sity gradients in the fuel. For this semi-analytic model 
of MagLIF, we assume temperature and density profiles 
that are comprised of two parts: a low-density hot spot 
region and a cold dense shelf region (see Fig. [J). This 
separation into two regions occurs when the preheat¬ 
ing is done rapidly across a region from r = 0 to Tpho, 
where Vpho < rg{tph)^ This separation occurs because 
the rapid preheating causes a blast wave to propagate ra¬ 
dially outward from rpho toward the outer radius of the 
overall fuel region, rg{t). This blast wave redistributes 
mass by boring out a low-density hot spot and pushing 
colder non-preheated fuel up against the liner’s inner sur¬ 
face. 

For the low density hot spot region, we assume a fuel 
(gas) temperature profile of the form 


Tg{r) = tAi 



ioi 0 < r < rh < r^, (105) 


and a radiation temperature given by 


Tr{r) = const. =Tb for 0 < r < r/i < Vg, (106) 


where Tc is the central (peak) temperature at r = 0, is 
the outer radius of the hot spot region (see Fig. H]), ^ is 
the power-law parameter that specifies the curvature of 
the profile, and Tb is the brightness temperature, which 
is found iteratively while calculating the radiative cooling 
rate, Pr (see discussion about Pr below, in Sec. imi. By 
comparing to detailed calculations (see for example Fig. 5 
in Ref. [l|) , we find the hot spot profile to be reasonably 
well modeled using ^ fv 2 with Nernst effects, and ^ ~ 
3.5 without Nernst effects. The results are not overly 
sensitive to this choice. We have found good agreement 
using anything from ^ = 2 to ^ = 8. 

Because of our isobaric assumption, specifying the tem¬ 
perature or the pressure at any point in the fuel automat¬ 
ically determines the other, i.e., Pg(r) ■ Tg{r) = const. 
Thus, Eq. 11051 means that the density profile in the hot 
spot region is given by 


Pgir) 




for 0 < r < r/j < Tg, (107) 


where Pc = pgTg/Tf. is the density at r = 0. 

In the cold dense shelf region, the fuel temperature and 
the radiation temperature are given by 

Tg{r) = Tr{r) = Tb ■ for < r < rg, (108) 


fJournal Reference: Phys. Plasmas 22 , 052708 (2015); http://dx.doi.Org/10.1063/l.4918953 
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and thus by our isobaric assumption, the fuel density in 
the shelf region is given by 

■ ( ~ ) for < r < Tg. (109) 
The temperature gradients in each region are thus given 

by 

^ ^ f-y (T’c-Tb) for 0 <r<r?,<rg 

for rh<r<rg. 

( 110 ) 

Note that the temperature and density profiles are nearly 
flat in the shelf region and that Tg{r), Tr{r), VTg, Pg{r), 
and Vpg are all finite throughout the entire fuel region. 
The reasons for using these functional forms, as well as 
how to find Tc, pc, rh, and Tb, are deferred to the dis¬ 
cussion below on radiative losses from the fuel fSec. lIUl) . 

With these profiles, we can now define the following 
radially dependent number densities in the fuel; 

ni = ni{r) = Ni-pg{r)/mg ( 111 ) 

Be = ne(r) = Zgm (112) 

Us = ns{r) = Ns ■ pg{r)/mg, (113) 

which are for the ions, electrons, and the ions of species s, 
respectively. Also, by assuming a “frozen-in” condition 
for the axial magnetic field in the fuel plasma, we have 
Bzg{r)/pg{r) = const. = Bzg/pg, and thus 

B^g = B^g{r) = ^Pg{r). (114) 

Pa 

Note that the Nernst effect (discussed below) actually 
breaks the frozen-in condition (c/. Fig. 5 in Ref. [l|). 
Nevertheless, this expression for Bzg{r) provides us with 
a simple first-order approximation of the radial depen¬ 
dence of the axial magnetic field that is accurate enough 
for the purposes of this model. 


J. Radiative losses 

The next term to calculate in Eq. [77] is the radiative 
cooling rate, Pr- This calculation determines the fuel 
temperature and density profiles described in Eqs. 1105b 
llOOI bv finding Tb iteratively while applying isobaric and 
conservation of mass arguments. The formulation is de¬ 
rived from a two-temperature gray model of nonequilib¬ 
rium radiative diffusion. Eor more details on this type of 
model, see the discussion on pages 261-262 in Ref. 

Using a two-temperature model was motivated after 
studying the results of full radiation magnetohydrody¬ 
namics simulations of MagLIE and observing the follow¬ 
ing: (1) the radiation temperature is nearly constant 
throughout the fuel, dropping only slightly in the shelf 
region as r approaches (2) the radiation temperature 


is significantly lower than the fuel temperature in the 
hot spot region; (3) the radiation temperature is roughly 
equal to the fuel temperature in the shelf region; (4) the 
radiative flux across the fuel-liner interface is well ap¬ 
proximated by (1 — as)aTs, where Ug > 0.9 is the albedo 
of the liner’s inner surface and Ts(r) is the temperature 
in the shelf region (both the fuel and radiation tempera¬ 
tures since they are equal in the shelf region). 

These phenomena occur because the liner’s inner sur¬ 
face heats and compresses rapidly to a point where it can 
no longer absorb or transmit radiation efficiently. Thus, 
this surface begins to reflect/reemit radiation back into 
the fuel region, effectively trapping a significant portion 
of the radiation in the fuel region. Moreover, thermal 
diffusion from the liner’s hot inner surface to the colder 
material deeper within the liner region is slow and on a 
timescale that is long compared to the implosion time. 
This slow diffusion is due in part to the high heat ca¬ 
pacity of the liner material, and it helps to enable the 
radiation trapping. 

These observations inspired the following two- 
temperature gray model for radiative losses over the vol¬ 
ume of the fuel, where opacity effects become stronger in 
and near the cold dense shelf region. Erom Eq. 6.62 in 
Ref. [tI, we have 

ig = AKraT^ - 4.Kp(TT^, (115) 

where Eg is the energy density of the fuel, is an av¬ 
eraged opacity described by Eq. 6.60 in Ref. [tI, and Kp 
is the Planck mean opacity. The first term on the right 
hand side of Eq. 11151 describes the rate at which the fuel 
absorbs energy from the radiation field, while the second 
term describes the rate at which the fuel loses energy to 
the radiation field. Since the radiation is trapped in the 
fuel region by the liner’s hot inner surface, and since, in 
the shelf region, the fuel mass and radiation fields are in 
thermal equilibrium, we assume that Kr ^ Kp, and thus 


ig = -AKpaTg 



(116) 


Note that a consequence of Eg. 11161 is that there is no net 
energy exchange between the fuel and the radiation field 
when the two are in thermodynamic equilibrium [T^ = 
Tg). Next, from Eq. 5.22 in Ref. It^. we have 


AKpaTg = J = Abr ■ Z'^gn^ne^/T'g, (117) 

where J is the frequency-integrated emission coefficient 
for the bremsstrahlung mechanism and Abr = 1-57 x 
10-40 ni3 . K-5 . J/s.^ Substituting Eq.|II7|into Eo.ITTbl 
and integrating over the volume of the fuel from r' = 0 
to r' = r gives the cumulative radiative cooling power 
from the fuel as a function of r 


Prv{r) = Abr-2TTh- Zg 



r'dr'. 
(118) 
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Note that since we have set = Tg in the shelf region, 
the contribution of the shelf region to the Prv integral 
is zero; this is consistent with the fuel’s mass and ra¬ 
diation field being in thermodynamic equilibrium in the 
shelf region. The first term in Eq. 11181 represents a radi¬ 
ation source and the second term represents a radiation 
sink due to opacity and absorption. Note that for the 
case of Tj. = 0 and constant density and temperature 
profiles, integrating Eq. 11181 over the entire fuel region 
gives the familiar bremsstrahlung radiation power from 
an optically-thin volume emitter 

Prv = Abr ■ ZgUifieyJ^g ■ ■ h. (119) 

Next, to find Tr^ we invoke the additional constraint 
that, in the shelf region (r^ < r < r^), the volume radi¬ 
ation described by Eq. 11181 must be consistent with the 
gray body surface radiation described by 

Prs{r) = (1 - Oi8)(yT^{r) ■ 2-Kr ■ h. (120) 

To do this, we first note that as long as Prv{rh) = 
Prsi^h), then Prv{r) = Prsi^) throughout the shelf re¬ 
gion because of the functional form assumed for the shelf 
temperature in Eg. 11081 i.e., Tg{r) = Trir) oc r~^P. This 
functional form was chosen to ensure that Eg. 11201 would 
return a constant value, and thus enforce our assumption 
of no radiative loss or gain within the shelf region. We 
also note that, in the hot spot region (0 < r < r?i), we 
have assumed Tr = const. = Tb, and thus Tr{rh) = Tb- 
Therefore, to ensure consistency in the shelf region, we 
must find the values of Vh and Tb that satisfy 

Prv{rh) = Prs{rh), ( 121 ) 

where 

/ Th 

'^i'^ey/Tg 

( 122 ) 

Prsirh) = - as)crTg ■ 2Trrh ■ h. (123) 

To find the values of and Tb that satisfy Eqs. 1121b 
mi we use a simple bisection method and evaluate 
the expressions for Prv(Th) and Prs(Th) iteratively un¬ 
til Prv{rh) = Prs{rh) to within 1%^ Also, before we 
can test our iterative guesses for Vh and Tb in Eos. 1121b 
11231 we need to find pc and Tc to complete the profile 
definitions given in Eqs. I105H1091 To do this, we use our 
isobaric assumption and conservation of mass arguments. 
Additional computational details for finding Tb, Pc, 
and Tc are provided in Appendix [B1 

Throughout a given simulation, Og = 0.9 is held fixed. 
However, the results are not very sensitive to this choice. 
In practice, we find that any reasonable value, say in the 
range of 0.5-0.95, works fine. 

Since, in the shelf region, Prvif) = Prsi^) = const. = 
Prvirh) = Prying), the total radiative cooling rate for 



Eq. [77] can be taken from anywhere in the shelf region; 
we take 

Pr=Prv{rg). (124) 

Note that in this model, if the entire fuel region is pre¬ 
heated uniformly, then the cold dense shelf region will 
never exist (only a hot spot region will exist), and thus 
all of the fuel mass will contribute to the radiation losses. 
Also, due to the isobaric assumption and the conserva¬ 
tion of fuel mass, the absence of a shelf region means 
that the central peak temperature, Tc, will be lower for 
the same amount of energy deposited in the fuel. For 
these reasons, uniformly preheating the entire fuel is not 
optimal. Moreover, preheating only a central portion of 
the fuel, say from r = 0 to Tpho ~ 0.5 • rg(tph), broadens 
the optimal MagLIF operating space by enabling the use 
of higher initial fuel densities (i.e., if these higher initial 
fuel densities were used in cases where all of the fuel was 
preheated, then the radiation losses would cool the pre¬ 
heated fuel too rapidly to obtain good fusion yield). This 
is fortunate since the experimental MagLIF program at 
Sandia is presently investigating the approach of preheat¬ 
ing the fuel via a laser with a small beam radius relative 

fO Cgitpfi). 

Finally, we comment that assuming this shelf region is 
present in experiments may be idealistic. Azimuthal or 
other 3D asymmetries (particularly during laser preheat¬ 
ing) might make it impossible to establish or maintain 
this shelf region. If a shelf region cannot be maintained 
(which would be difficult to diagnose experimentally), 
then the optimal operating space for MagLIF might not 
be as broad as that suggested by simulations where the 
initial fuel density is scanned over and rpho < rg{tph). 


K. Magnetized electron and ion thermal conduction losses 


The next term to calculate in Eq. |73is the energy loss 
rate due to electron and ion thermal conduction. To do 
this, we use the Epperlein-Haines transport equations for 
the electrons^! and the Braginskii transport equations for 
the ionsi^ The energy loss rate due to electron thermal 
conduction is given as a function of r by 

dT 

Pce(r) = 2'Krh ■ tle{Xe) * , (125) 

or 

where 

^ _ UekTgTei _ 6.18 -b 4.66Xe _ ,^20') 

® me 1.93-b 2.31*6 + 5.35a;e-b Xg 

is the coefficient for electron thermal conduction perpen¬ 
dicular to B^gj *e = WceTei is the electron Hall parameter, 
<^ce = QePzgl'iTT'e is the electron cyclotron frequency, and 
Tei = 1 /vei is the average time between electron-ion col¬ 
lisions. Here, Vei is the electron-ion collision frequency 
given by^ 


4v^ {Y^.nsZl) g^lnA 
{47reof 3v^(fcTg)'/^ 
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where we have accounted for multiple ion species, s, and 
where the In A calculations are described in Appendix El 
Similar to electron thermal conduction, the energy loss 
rate due to ion thermal conduction is given as a function 
of r by 


where 


Pci{r) = 

dT 

2Krh • ni{xi) • 

or 

(128) 

TlikTgTii 

2.645 -b 2x1 

(129) 

mi 

0.677-b 2.70a;2 -b 4 


is the coefficient for ion thermal conduction perpendic¬ 
ular to Bzg, Xi = uJciTii is the ion Hall parameter, 
Wei = geZgBzglfhi is the ion cyclotron frequency, and 
Tii = is the average time between ion collisions. 

Here, va is the ion-ion collision frequency given by^ 


_ 40Fni,Z^gf InA ( n,Zl \ 

{47reo}^ 3v^ [kTgf^^ nhZl j ’ 

(130) 

where Uh = Ud + nt, Z^ = Zd = Zt, and where we 
have assumed that the plasma is comprised mainly of fuel 
deuterons and tritons (subscript h for hydrogen isotopes), 
with an average dominant ion mass of 


rrih = 


Nduid + Ntmt 
Nd + Nt 


(131) 


that is much less than the mass ms^d,t of any of the 
dopant/contaminant particles. The In A calculations are 
described in Appendix El 

Finally, the total radially dependent energy loss rate 
due to thermal conduction is given by 


Pcir) = Pce{r)+P„ir). (132) 


For this calculation, we use radially dependent param¬ 
eters ni(r), rieir), Xe(r), and Xi{r) [and thus Bzg(r), 
Vei{r), Vii(r), nh{r), ns(r), and InA(r)], since they are 
all readily calculable from our given expressions for Tg(r) 
and Pg{r). 

While studying thermal conduction using full radiation 
magnetohydrodynamics simulations (as well as our hot 
spot model), we observed that ion thermal conduction 
dominates over electron thermal conduction in the inner 
regions of the fuel hot spot. By contrast, near the edge of 
the hot spot (i.e., r k, r^), the thermal transport is dom¬ 
inated by electron conduction. Thus there is a handoff 
from ions to electrons at some radius in the fuel hot spot 
region. Because of this, and because our hot spot pro¬ 
file is prescribed and determined by the radiation model, 
we wish to be conservative with regards to the expected 
losses due to thermal conduction, and therefore we look 
for a maximum in Pc{t). That is, for the thermal con¬ 
duction loss power from the hot spot region to the shelf 
or liner region, we use 

Pch = max |A;(?’)| for 0 < r < r?, < rg. (133) 


For the thermal conduction loss power from the shelf re¬ 
gion to the liner region, we simply use 

Pcs = Pcirg). (134) 

The overall thermal conduction loss power from the fuel 
to the liner depends on whether or not the shelf region 
exists, and thus, for the thermal conduction power in 
Eq.jTTl we use 


Pc = 



for rn < Vg 
for rh=rg. 


(135) 


Note that even if a shelf region exists, we still need to 
calculate Pch in order to estimate the shelf erosion rate 
(see Eq. 11511 below). 


L. End losses 

The final term to calculate in Eq. |77]is the energy loss 
rate due to fuel losses out of the ends of the imploding 
cylinder. Bends- We estimate this loss by assuming that 
the fuel (with its associated mass and energy densities) 
flows across the top (z = h) and bottom (z = 0) planes of 
the imploding region at the hydrodynamic sound speed, 
which is determined by the pressure and mass density of 
the fuel. Therefore, the energy losses out of the top and 
bottom planes, respectively, are given by 

. P rctopii) 

Etop = (3/4)^ * TT / ■ dr 

^9 Jo 

Ebot = (3/4)"^ • / Cg{r)- 2'Kr ■ dr, 

Jo 

where rtop(t) and rbotp) are the radii for the top and 
bottom apertures that the fuel can escape through, and 
Cg(r) is the hydrodynamic sound speed, given by 

Cg{r) = \JlgPg/Pg{r), (138) 

where 7 g = 5/3 is the ratio of specific heats for an ideal 
gas. The factor of (3/4)'^ is a result of the derivation 
presented in Ref. [l|. Finally, the total energy loss due to 
end losses is2^ 


(136) 

(137) 


Bends — Etop T Ebot- (139) 

Similarly, for mass losses from the fuel, we assume 
ion density flows across the apertures, and thus, for ion 
species s, we have 


j-Ctopii) 

Ns, top = (3/4)^ / ns{r)cg{r) ■ 27Tr ■ dr 
Jo 

l-rhoid) 

Ns, hot = (3/4)^ / ns{r)cg[r) • 2'Kr • dr 

Jo 

Es,ends — E^s,top T E^s,bot- 


(140) 

(141) 

(142) 
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The number of deuterons and tritons in the fuel change 
because of fusion reactions as well as end losses (see 
Eos. 11631 and ITMl belowl. By contrast, end losses are 
the only way for prescribed dopant and/or contaminant 
particle numbers to change. Thus, we have 

^mix — E Ns ,ends 1 (143) 

s^d,t 

where Nmix{t) = t N^s{t) is the total number of 

dopant and contaminant particles. Note that we only 
need one differential equation to describe the evolution 
of all of the “mix” particles because their number dis¬ 
tribution amongst themselves remains constant in time, 
i.e., 

fs^d,t{t) = = const. = fs{to), (144) 

^^mix } 

and therefore 

Ns^ddit) = (145) 

Finally, should fuel preheating be done with a laser 
entering from above the imploding region, then these 
end loss calculations should be made for at least the top 
plane, where the radius of the top aperture is given by 

rtop{t) = min {tlec, (146) 

and where tlec > Ci, is the radius of the laser en¬ 
trance channel. The calculation for the bottom plane, 
however, depends on whether the bottom plane is com¬ 
pletely closed off with electrode material (i.e., an implo¬ 
sion “glide plane”) or if it has an opening for a beam 
dump. Should a beam dump be used, then the radius of 
the bottom aperture is given by 

rbot{t) = min{rdu„p,rg(t)}, (147) 

where rdump is the radius of the beam dump. 


Note that, like thermal conduction, Nernst losses also de¬ 
pend on a function of the Hall parameter, F{Xf. = uiceTei), 
as well as gradients in the fuel temperature. Also, we 
assume that the magnetic flux lost by the fuel is trans¬ 
ported to the liner, and thus 

= <^zlO + [<^zgO-<^zg{t)]. (150) 

As the axial flux is transported into the liner, it is likely 
that it is dissipated rapidly by resistive diffusion, which 
would heat the liner material somewhatHowever, the 
energy density in the liner for either case (either pure 
magnetic energy or ohmic dissipation leading to thermal 
energy) is roughly the same, and thus will affect the liner 
dynamics similarly given the other simplifying assump¬ 
tions made in our model. Thus, for simplicity, we ignore 
ohmic dissipation of the axial magnetic field. Finally, we 
assume that the axial magnetic flux in the vacuum region 
is conserved, thus ^zv{t) = ^zvo- 


N. Erosion of the fuel’s dense outer shelf region 

The cold dense shelf region, generated by preheating 
only a central portion of the fuel, provides a buffer region 
in the fuel between the hot spot and the cold liner wall. 
This buffer region significantly reduces radiation losses 
by effectively reducing the amount of fuel mass that con¬ 
tributes to the radiation losses. It also reduces thermal 
conduction losses and axial magnetic flux losses from the 
fuel to the liner because the temperature gradient at the 
edge of the nearly flat shelf region is much less than that 
at the edge of the hot spot region. However, this buffer 
region is not static. It begins to erode away immediately 
after its formation due to thermal conduction from the 
hot spot region to the shelf region. To describe this ero¬ 
sion, we differentiate E = ^NkT with respect to time, 
rearrange the terms, and make simple substitutions to 
approximate the mass transfer rate from the shelf region 
to the hot spot region as 


M. Magnetic flux loss due to the Nernst thermoelectric 
effect 


. _ 2 nii {Ecb Pcs) 

" 5{1 + Zg) kiTh-TsY 


(151) 


With all of the terms in Eq. [77] thus defined, the next 
important quantity to calculate is the flux loss of the axial 
magnetic field from the fuel region to the liner region due 
to the Nernst thermoelectric effect. For MagLIF with 
a hot fuel region, flux losses due to the Nernst effect 
dominate over flux losses due to resistive diffusion alone 
To model flux losses due to the Nernst effect, we follow 
Ref. m to write 




T{Xe) 


—2'Kr ■ T (xe) 


k dTg 
Qe dr 


1.5Xg -I- 3.053a;e 
xt -f 14.79a;2 -h 3.7703’ 


(148) 

(149) 


where Th is the average temperature in the hot spot re¬ 
gion and Ts is the average temperature in the shelf re¬ 
gion. This states that thermal energy must be deposited 
in the shelf for a particle to make a “quantum” jump 
from having an average thermal energy of ^kTs in the 
shelf region to an average thermal energy of | fcT/ in the 
hot spot region; the larger the temperature jump and/or 
mass transfer rate, the larger Pch — Pcs needs to be. The 
average temperatures in each region are found by invok¬ 
ing our isobaric assumption, i.e., 

fh = Pgfg/ph 

Ts = PgTgj Ps, 
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where 


Ph = fmhrrig/ {irrlh) (154) 

Ps = (1 - fmh) rUg/ [tt {rl - rl) h] , (155) 

are the mean fuel densities in the hot spot and shelf re¬ 
gions, respectively, and fmh is the fraction of the total 
fuel mass that is part of the hot spot. We keep track of 
fmh rather than the absolute fuel mass in the hot spot 
since the absolute mass changes (due to fusion reactions 
and end losses). Thus we have 


fmh = ms^hlmg. (156) 

Note that initially (at the time of preheat), we have 


fmhi^iph') — 


'^phO 

'g{tph) 


1 2 


(157) 


where if laser absorption is being calculated, Vpho = rb. 
From our example above using the laser preheating pa¬ 
rameters of recent MagLIF experiments^ we get 


fmh{tph') 


0.492 

2.325 


[0.21]^ = 0.04. 


(158) 


Thus, initially, only about 4% of the fuel mass is part of 
the hot spot region, though, radially, the hot spot region 
occupies 21 % of the overall fuel radius Vg, and quickly 
expands to occupy ~ 90% of rg before the shelf begins 
eroding away with any signihcance. 


O. Primary fusion reaction rates 


TABLE III. Coefficients for Eos. 11621 (from Refs. and 1^1 . 



DT 

DD,®He 

DD,T 

Co 

6.6610 

6.2696 

6.2696 

Cl X 10^2 

643.41 X 0.98^ 

3.5741^ 

3.7212 

C2 X 10^ 

15.136 

5.8577 

3.4127 

Cs X 10® 

75.189 

7.6822 

1.9917 

C 4 X 10® 

4.6064 

0 

0 

C 5 X 10 ® 

13.500 

-0.002964_^ 

0.010506 

Ce X 10® 

-0.10675 

0 

0 

C 7 X 10® 

0.01366 

0 

0 


^ For DT, we have multiplied the Ci coefficient of Ref. by 0.98 
for a slightly better fit to the peak of the tabulated {av) data 
published in Ref. [93. 

For DD,®He, and for T > 100 keV, we set C 5 = 0 to avoid 
{<yv)dd 3He Too near 965 keV. Though unlikely for any 
reasonable MagLIF case, this ±00 has caused computational 
problems when scanning over a large parameter space, where 
very low density plasmas can be heated to 1 MeV. This 
replacement causes a discontinuity at 100 keV; to remove this 
discontinuity, we also set Ci X 10^^ = 3.5741 X 1.0172 for 
T > 100 keV. 


where the coefficients Cq-C^ are the fitting parameters 
provided in Table uni 

The reactivity rates, along with end losses, provide the 
rate of change of the total number of deuterons and tri¬ 
tons in the fuel, respectively, as 


Nd = —Ndt — ‘2-NddpHe — “^Nddd — Nd^ends ( 163 ) 
Nt = —Ndt — Nt^ends- (164) 


For an arbitrary deuterium to tritium fuel ratio, the 
DT reaction rate is given by^i 

Ndt = h ndnt{<Tv)dt ■ ■ dr, (159) 

Jo 

and the two dominant DD reaction rates are given by^i 
h po 

NddPRe =2 nl{av)ddPRe ■ Sirr • dr (160) 

Ndd,t = 2 J '^d{^'^)dd,t ■ 27rr • dr, (161) 

where rid and ut are the temporally and radially depen¬ 
dent deuteron and triton number densities in the fuel, as 
defined by Eq. 11131 with s = d and s = t, respectively, 
and where {av)dt, {crv)dd,me, and {av)dd,t are the tem¬ 
porally and radially dependent reactivity parameters for 
the DT, DD-^He, and DD-T reactions, respectively. The 
reactivity parameters are calculated usin g^^i^^ 

{av) = exp (-3C^/^C) (162a) 

C2Tg,heV + 

1 -I- CsTg^keV + 

(162b) 

i = Co/Tl'lv^ (162c) 


The reactivity rates also provide the total fusion power, 
as given by 

Pf — NdtQdt T NddpHeQddpHe T Ndd,tQdd,t-) (165) 
where the energy yields per DT, DD-^He, and DD-T re¬ 


action are, respectively, 

Qdt = 17.6 X 10® • Qe (166a) 

Qdd,me = 3.27 X 10® • (166b) 

Qdd,t = 4.03 X 10® • ge- (166c) 

Since one neutron is released per DT reaction and one 
neutron is released per DD-^He reaction, the primary 
DT and primary DD neutron yields are given by 

Ydt,n = Ndt (167) 

Ydd,n = Ndd,mei (168) 

and thus the total primary neutron yield is 

Yn = Ydt,n + Ydd,n- (169) 

Finally, the total fusion energy yield is 

Y = NdtQdt + Ndd,meQdd,me + Ndd,tQdd,t- (170) 
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P. Summary of the semi-analytic MagLIF model 


MODEL VERIFICATION 


The dynamics of our model are described by Eqs. Emu 
[Ml lik [371 [ 5 ^ [771 [T4^ [T4^ fTKm [T591 [TO fWl [T^ 
and 11641 These 2Nis + 13 ordinary differential equations 
are repeated here for emphasis and clarity: 


^oc 

L 


(171a) 

^c/ Rloss 

c 


(171b) 

Pc — Ii{Lv P Lie) 


(171c) 

Lq P Ly P Lie 


Pl,s=i — 1 Pl,s—i 

27rn,i • h 

(171d) 

mis 


(* = 2,3,.. 

.,Nu-l) 

PgPPB,^ - Pi,S=1 

misl2 

• 2'Krg ■ h 

(171e) 


Pl,s=Nu - PBeiv - Pb,^ o , 
ri = -—-27rr/ • h 

TOis/2 


(171f) 

(171g) 


+ [Pr + Pc + Pq. — Pbb) /Nls 

is = l,2,...,Ni,) 


Pg — PpdV P Pph P Pa Pr Pc Pends 


Nrr 




— ^ ^ ^s,ends 

s^d,t 

Qe or J 


fmh — / BXg 


Ndt = h / ndnt{(jv)dt ■ ^irr ■ dr 
Jo 

h pH 

n I 2 


P^ddpHe — 2 


Nf = -Ndt - N, 


(171h) 

(171i) 

(171j) 

(171k) 

(1711) 

(171m) 


nd{crv)ddPHe ■ 27rr • dr 

h 

Ndd,t = o / '^d{<^'^)dd,t ■ 27rr • dr (171n) 

2 Jo 

Nd = —Ndt — 2Ndd,^He — ‘JNdd,t — Nd,ends (171o) 


ends' 


(171p) 


Note that if the liner current is prescribed, rather 
than derived from the voltage-driven circuit model, then 
Eqs. ll71al[T7Tcl are discarded, and only 2Nis -b 10 ordi¬ 
nary differential equations are required to describe the 
dynamics of our model. Furthermore, if we ignore the 
Nernst effect, preheat all of the fuel uniformly, use ei¬ 
ther pure D 2 fuel or a 50/50% DT mix, and we ignore 
dopants and contaminants, then the number of required 
differential equations drops to only 2Nis + 5. 


In Fig. El we present SAMM simulation results where 
the intent was to reproduce the ID results presented in 
the original MagLIF paper These SAMM results show 
that our simplified model does indeed capture the ID 
behavior presented throughout Ref. [l|. All of the simula¬ 
tions presented in Fig. E] used DT fuel. 

In Fig. EKa), we have plotted the fuel radius, liner ra¬ 
dius, drive current, and average fuel temperature as a 
function of time for the preliminary point design dis¬ 
cussed in Ref.Jh. This figure should be compared with 
Fig. 4 in Ref. [J The SAMM simulation produced a fu¬ 
sion energy yield of 970 kJ, compared to about 500 kJ in 
Ref. [l|, and resulted in a maximum convergence ratio of 
25, which matches the preliminary point design of Ref. [l|. 

In Fig. Elb), we have plotted the normalized fuel tem¬ 
perature, magnetic field strength, and fuel density at 
stagnation (at peak burn rate) as a function of the nor¬ 
malized fuel radius for two SAMM simulations, one with 
the Nernst effect included (solid lines), and one without 
the Nernst effect included (dashed lines). All quanti¬ 
ties are normalized to the simulation without the Nernst 
effect 1 ^ This figure should be compared with Fig. 5 in 
Ref. [l|. The yield for the simulation with Nernst is 72% 
lower than the simulation without Nernst (compared to 
about 70% in Ref. [l|). The simulation with Nernst re¬ 
sulted in an axial magnetic flux loss of about 64% (com¬ 
pared to about 70% in Ref. [l|). 

Figures EKc) and EKd) should be compared with 
Figs. 6(a) and 6(b) in Ref. [l|. In Ref. [l|, the initial fuel 
densities that result in solutions with convergence ratios 
of 10, 20, and 30, were found as a function of the initial 
fuel temperature. Using these initial conditions directly 
with SAMM, shown here in Fig. EKd), we generated the 
results shown in Fig.E^c). The resulting convergence ra¬ 
tios were 16.3±1.4, 25.7±3.4, and 32±5.5, and thus are 
only approximately equal to 10, 20, and 30. 

In Fig. EKe), we have plotted the optimum fusion en¬ 
ergy yields per unit liner length as a function of the ini¬ 
tial axial magnetic field strength. Here we have held the 
maximum convergence ratio fixed at either 10, 20, or 30. 
To find these optimum constant convergence ratio solu¬ 
tions, we had to scan over the input parameter space that 
consists of the initial preheat energy and the initial fuel 
densi^. These results should be compared with Fig. 7 in 
Ref. [3 

Figures EDO and EKg) should be compared with 
Figs. 8(a) and 8(b) in Ref.[ll In Fig.EJf), we have plotted 
the optimum fusion energy yields per unit liner length as 
a function of the peak drive current for two cases, one 
with a-deposition turned on, and one with a-deposition 
turned off. Here, as in Ref.[l|, for the case with a-heating, 
we scanned over the input parameter space that consists 
of the initial preheat energy and the initial fuel density to 
find the optimum solutions that had a maximum conver¬ 
gence ratio of 20 (for the case with no a-heating, the same 
input parameters were used, but a-deposition was simply 


1 Journal Reference: Phys. Plasmas 22, 052708 (2015); http://dx.doi.Org/10.1063/l.4918953 
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(a) cf. Fig. 4 in Ref. 1 


(e) cf. Fig. 7 in Ref. 1 


(i) cf Fig. 11 in Ref. 1 
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(d) cf Fig. 6(b) in Ref. 1 



10 


o 10 


■a 
1 10 
>- 
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FIG. 4. Simulation results from our semi-analytic M^LIF model (SAMM), where the intent was to reproduce the simulation 
results presented in the corresponding figures of Ref. [J. 
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turned off, thus the convergence ratios varied somewhat). 
In Fig. EDg), we have plotted the ratio of the peak fuel 
temperature with a-heating turned on to the peak fuel 
temperature with a-heating turned off. We have also 
plotted the ratio of the total fusion energy yield to the 
total energy absorbed by the target (liner and fuel), i.e., 
the “target gain”. 

In Fig. |4l)h), we have plotted the optimum fusion en¬ 
ergy yields per unit liner length as a function of the initial 
liner aspect ratio (ARl. This figure should be compared 
with Fig. 10 in Ref. [l[ In the SAMM results, the tran¬ 
sition to higher yields with higher aspect ratios is a bit 
more abrupt, and then rolls off sooner; however jthe gen¬ 
eral trend is the same as that shown in Ref. [J. Also, 
although not explicitly shown in Ref. [l|, the eventual 
rolloff in yield with higher aspect ratio is expected. That 
is, to keep the implosion time constant with higher as¬ 
pect ratio, the overall liner mass must be decreased while 
the liner’s initial outer radius is increased (and while the 
liner’s initial wall thickness is decreased). Eventually, the 
liner’s mass decreases to a point where it cannot provide 
sufficient tamping and inertial confinement at stagnation. 

In Fig. EKi), we have plotted the normalized yield 
as a function of the premixed fraction of beryllium in 
the fuel. The yields are normalized to the case with 
no mix/dopants. This figure should be compared with 
Fig. 11 in Ref. [l|. 

In Fig.EKj), we consider laser propagation through the 
fuel using the analytic treatment discussed in Sec. IIIGI 
Here we plot the laser heating front propagation distance 
as a function of time for the preliminary point design of 
Ref. [l|. The laser pulse duration is 10 ns. At 5 ns (red), 
the front of the “bleaching wave” has reached 0.56 cm 
into the fuel, which is consistent with Fig. 13(b) of Ref. [l|. 
By 8 ns (blue) and 11 ns (green), the bleaching wave has 
propagated completely through the 0.65-cm-long implod¬ 
ing fuel region, which is consistent with Figs. 13(c) and 
13(d) of Ref. [H, respectively. 

In Fig. EKk), we present a color-filled contour plot of 
total fusion energy yield as a function of both the fuel 
preheat temperature and the laser beam spot size (ra¬ 
dius). Overlaid on top of the color-filled contour plot 
are the contours of constant preheat energy required to 
obtain the preheat temperature at the given spot sizes. 
This fi^re should be compared with Fig. 15 in Ref. [l|. 
In Ref. U, this figure was generated by additionally scan¬ 
ning over a range of initial fuel densities to obtain only 
the solutions where the convergence ratios were 25. For 
simplicity, we used the input densities and preheat ener¬ 
gies found by this previous effort to generate the SAMM 
results presented here in Fig. EKk), and thus the conver¬ 
gence ratios were allowed to vary somewhat. 

Finally, in Fig. EKl), we have plotted results for mass 
end losses out of the fuel region through the laser en¬ 
trance hole using the analytic treatment discussed in 
Sec.imi Here we are plotting the ratio of mass remain¬ 
ing after fusion burn to the mass of the fuel at the start 
of the simulation as a function of the axial length of the 


liner. We have done this for three cases, where the ra¬ 
tio of the radius of the laser entrance hole to the initial 
radius of the fuel is either 1.0, 0.5, or 0.25. This figure 
should be compared with Fig. 16 in Ref. [H 


IV. SUMMARY, FUTURE WORK. AND CONCLUSIONS 

In this article, we have presented a new semi-analytic 
model of MagLIF. This model, implemented in a code 
called SAMM, accounts for: (I) preheat of the fuel (op¬ 
tionally via laser absorption); (2) pulsed-power-driven 
liner implosion; (3) liner compressibility with an ana¬ 
lytic equation of state, artificial viscosity, internal mag¬ 
netic pressure, and ohmic heating; (4) adiabatic com¬ 
pression and heating of the fuel; (5) radiative losses and 
fuel opacity; (6) magnetic flux compression with Nernst 
thermoelectric losses; (7) magnetized electron and ion 
thermal conduction losses; (8) end losses; (9) enhanced 
losses due to prescribed dopant concentrations and con¬ 
taminant mix; (10) deuterium-deuterium and deuterium- 
tritium primary fusion reactions for arbitrary deuterium 
to tritium fuel ratios; and (11) magnetized a-particle 
heating. We have also shown that SAMM is capable of 
reproducing the general ID behavior presented through¬ 
out the original MagLIF paper (i.e., Ref. [l|). 

Presently, we are using SAMM to further explore the 
parameter space of MagLIF. These efforts are comple¬ 
mentary to other studies presently taking place using 
full ID, 2D, and 3D radiation magnetohydrodynamics 
codesj ^^'^^'^^ In particular, we are using SAMM to study 
the parameter space surrounding the experimental plat¬ 
form of Ref. [l^ (i.e.. an initial axial magnetic field of 10 T, 
a preheat energy of less than 2 kJ, an initial fuel density 
of about I mg/cm^, and a peak drive current of about 
20 MA). We are also studying how this parameter space 
could change as a series of planned upgrades are made 
to the Z facility over the next few years. These upgrades 
are intended to bring MagLIF’s experimental platform 
closer to the parameters of the preliminary point design 
described in Ref. [l| (i.e., an initial axial magnetic field of 
30 T, a preheat energy of 8 kJ, an initial fuel density of 
3 mg/cm^, and a peak drive current of 27 MA). Addi¬ 
tionally, we are studying how this parameter space could 
change as the peak drive current is increased to about 
60 MA, in consideration of MagLIF on possible future 
pulsed-power accelerators, such as the conceptual Z-300 
and Z-800 acceleratorsThe results of these parameter 
space studies using SAMM will be presented in a future 
publication. 

In closing, we note that the development of SAMM has 
led to several physical insights that were not fully appre¬ 
ciated previously (e.g., the dependence of radiative loss 
rates on the radial fraction of the fuel that is preheated). 
Additionally, this model’s accessible physics and fast run 
times (^ 30 seconds/simulation) make it a useful peda¬ 
gogical tool, especially for students, experimentalists, or 
any researcher interested in MagLIF. 


1 Journal Reference: Phys. Plasmas 22, 052708 (2015); http://dx.doi.Org/10.1063/l.4918953 
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Appendix A: Calculations of the Coulomb logarithm, In A 


Following Ref. l80|, the Coulomb logarithm is given by 

Ad 


In A = max < 1, In 


(Al) 


^qva ^ --- 

2meUT, 


(AS) 


is the quantum mechanical minimum impact parameter. 
Here, h = 1.05457 x 10“^^ J • s is the reduced Planck 
constant and 


VT^ = ^j2kT/me (A6) 

is the electron thermal velocity. 

For the In A calculation in Eq. |SS1 we use the following 
assignments: Zs Znuc,s, kT -)> 13.6 • qe ■ J2s^nuc,s^ 

Ue fie, and Us ^ fls. 

For the In A calculation in Eq. 11011 we use the following 
assignments: kT —>■ kTg, rie —?> he, and Us ^ fis- 

For the In A calculations in Eqs. 11271 and 11301 we use 
the following assignments: kT —^ kTg(r), Ug —>■ ne{r), 
and Tis risir). 


Appendix B: Additional computational details for finding 

Th, Tb, Pc, and Tc 

1. When the shelf region is not present 

This case occurs when either all of the fuel is pre¬ 
heated uniformly, or when the shelf region has completed 
eroded away (see Sec. IllNp . For this case, Vh = Vg, 
and we iteratively guess for Tb in our bisection algo¬ 
rithm. After evaluating the expressions in Eqs. 1121111231 
if PrsiTh) < Prvifh), then Tb is too low, and the cor¬ 
rect value for Tb will be found between this value that is 
too low and the average temperature of the entire fuel, 
Tg. [Note that we take Tg as the initial upper limit for 
Tb in the bisection algorithm because with rt = Vg and 
Tb = Tg, the resulting fuel temperature and density pro¬ 
files are flat, and Pryirh) = 0 while Prs{rh) > 0.] Con¬ 
versely, if Prs{rh) > Prvifh), then Tg is too high, and 
the correct value for Tb will be found between this value 
that is too high and zero. [Note that Tb = 0 results in 
Prs{rh) = 0 while Prv{rh) > 0, thus we take zero as the 
initial lower limit for Tb in the bisection algorithm.] 


where Ad is the Debye length of the plasma and bmin 
is the minimum impact parameter of the plasma. The 
Debye length is given by 


Xd = 


eokT 


{ne + Es 


(A2) 


The minimum impact parameter is conventionally taken 
as 


where 


2. When the shelf region is present 

For this case, we must ensure that rh and Tb are chosen 
such the following mass conservation argument is met: 


‘•gs 


= {i^-fmh)mg= / Pg{r) ■2-Kr ■ h-dr (Bl) 
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Here fmh is the fraction of the total fuel mass that is 
in the hot spot region, which is the system variable de¬ 
fined by Eqs. I156H157I (see also Ea. ll71kll . Also, we have 
used Eg. 11091 for Pg{r) and, in the last equality, we have 
used our isobaric assumption that PcTc = PgTg. Equa¬ 
tion IB3I provides the analytic relationship between Tb 
and Th that ensures mass conservation in the shelf region. 
However, if we were to iteratively guess for Tb in our bi¬ 
section algorithm, then solving Eq. IB3I for Vh in terms of 
Tb would require a numerical treatment. Therefore, we 
instead guess iteratively for r^, which must reside in the 
domain 0 < < rg (thus, for our bisection algorithm, we 

use r = 0 and r = Vg as our initial lower and upper lim¬ 
its, respectively). For each guess. Eg. IB3I provides the 
unique value of Tg that will ensure mass conservation. 
Note that TB{rh) is a monotonically decreasing function 
(i.e., Tb —>■ oo as —>■ 0, and Tb —>■ 0 as —>■ Vg). 
Thus, in our bisection algorithm, if Prs{rh) < Prvirh), 
then Tb is too low and Vh is too large; the correct value 
for Th will be found between this value that is too large 
and zero. Conversely, if Prsifk) > Pryi^h), then Tb is 
too high and is too small; the correct value for will 
be found between this value that is too small and Vg. 


3. Updating pc and Tc 


Before we can test our iterative guesses for Vh and Tb 
in Eqs. I121H1231 we first need to find pc and Tc to com¬ 
plete the prohle definitions given in Eqs. I105H1091 To 
do this, we use our isobaric assumption to first update 
Pg{i^h) = PgTg/TB = Pb- Next, we find the value for 
Pc that makes the total mass in the hot spot, as de¬ 
fined by Eq. 11071 consistent with mgh = fmh'm'g (see 
Appendix IB 41 for additional details). Then, we again ap¬ 
ply our isobaric assumption to get Tc = pgTgjpc- Finally, 
with all profile parameters defined and updated, we eval¬ 
uate Prvirh) and Prs{rh) and test whether the results 
meet our desired solution accuracy. 


4. Additional details for updating pc. 

The central mass density, pc = Pg(r = 0), is updated 
after updating Tb, r^, and Pg{rh) = Pb- The updated 
Pc value is found by ensuring that the following mass 
conservation argument is met: 

j-rh 

mgh = fmhmg= / pg{r) ■ 2nr ■ h ■ dr, (B4) 
Jo 

where Pg{r) is given by Eq. 11071 If the expressions given 
for Tg{r) and Pg{r) in Eos. ll05l and ll07l were permitted to 
extend to radii beyond r^, then we would have Tg(r) — >■ 0 
and Pg{r) —>■ oo when 


r 


rh 



(B5) 


Using r* to define the dimensionless radius x = rfrt., the 
hot spot mass density given by Eq. 11071 becomes 

Pg{r) = Pc (1 - x^) ^ ■ (B6) 

Evaluating the dimensionless radius x at r = Vh gives 

Xh=rhlr^, (B7) 


and thus the total fuel mass in the hot spot region, as 
given by the right hand side of Eq. IB41 becomes 


rugh = 



Pc 

(1 - a:«) 


• • 2'kx ■ h ■ dx 


= Pc ■ -xrl 



2x 

(l-a;«) 


• dx 


(B8) 

(B9) 


= Pc 


— \-h- 


2x 


= Pc ■ 'xrlh ■ 


(1-x?) 
2x 


■ dx 


(l-x«) 


dx. 


Since the average mass density in the hot spot 
Ph = ’nigh/i'xr'lh), we have 


(BIO) 

(Bll) 
region is 


Ph — Pc 


2x 


Now, from Eqs. | 


land 


{1-xi) 
we have 


• dx. 



(B12) 


(B13) 


where we have again applied our isobaric assumption; in 
this case, p^c = PbTb ^ Tb/Tc = pdP b- Solving for 
Pc in Eq. IB13I gives 


Pc = PB (l - 4) ■ (B14) 

Substituting Eo. IB 141 into Eo. lB12l and solving for Ph/P b 
gives 


Eh. 

Pb 


1-4 r 

Jo 


2x 
1 — 


• dx = V{xh), 


(B15) 


where Pixh) is a dimensionless, monotonically decreasing 
function of the dimensionless variable Xh', it is defined 
over our physically relevant domain of 0 < < 1 with 

a range of 0 < V{xh) < 1 (see Fig. [5]). Because P{xh) 
is a dimensionless function of a dimensionless argument, 
we only need to evaluate it once and store the results 
for continuous use throughout the simulation. With this 
done, we take the inverse of P(xh), with input argument 
Ph/PB, to get 


(B16) 

Finally, from Eq. IB14I we get the central density, pc, 
repeated here for completeness as 

Pc = PB (l - 4) ■ (B17) 
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FIG. 5. Plot of V{xh)^ the dimensionless, monotonically de¬ 
creasing function of the dimensionless variable as defined 
by Eq. [BT51 
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